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Lattice diffusion and internal local magnetic field gradients in solids are investigated by numeri- 
cal simulation of nuclear mangetic resonance(NMR) spin echo experiments. The Fourier-spectrum 
method is employed in order to solve the Bloch-Torrey equation with arbitrary magnetic field gra- 
dients in one- and two-dimensional lattice restrictions. 
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I. INTRODUCTION 

Transport phenomena are often encountered in solu- 
tions, fluids, gels, and quasi particles associated with 
elementary excitations in solids, and so on. In par- 
ticular, low dimensional transport phenomena are re- 
ceived much attraction in various research fields of elec- 
tronic devices such as polymeric conducting materials 1 , 
quasi one-dimensional organic conductors- and recharge- 
able batteries^. In these matarials or devices, low di- 
mensional diffusive or hopping motion such as charged 
ions, defects, quasi-particles(e.g., exciton, polaron, con- 
formon) in their crystalline lattices plays an important 
role in origins of their physical properties. Some of such 
molecular motions in crystalline solids are called "lat- 
tice diffusion," which is defined as transport phenomena 
of (quasi-)particles of which positions are restricted in 
lattice points in crystals. The motion often leads to mo- 
tional averaging of local magnetic fields in solids. Nuclear 
magnetic resonance(NMR) spectroscopy is one of promis- 
ing nondestructive methods for investigating these phe- 
nomena. Several attempts have recently been reported 
to make NMR measurements of ultra slow diffusion in 
solids&&±2. 

Among these, the pulse field gradient(PFG) technique 
with strong magnetic field gradient (MFG) has been given 
much attention for detection of such diffusions. How- 
ever, the PFG technique involves some problems in cases, 
which are (i) attenuation of sensitivity in association of 
application of PFG, (ii) a lack of magnitude of PFG for 
ultra slow diffusion measurements, and (iii) the case that 
the system has large internal MFG(G;) compared to ex- 
ternal PFG and/or large heterogeneous MFG. 

The problem (i) and (ii) are often pronounced in molec- 
ular diffusion measurements. Since the problem (iii) is 
not problematic in cases of systems with large diffusion 
coefficient, it has often been ignored. In such cases, the 
internal MGF is thought to be averaged to zero, and then 
effects of the internal MFG will be small or negligible. 
However, in solids, the diffusion coefficient can be quite 
small and a local internal MFG, which in some cases het- 
erogeneous and/or anisotropic, can be remained^. 

In such situations, an effective MFG is no longer 
the same as the external MFG (G c ), which is likely to 



mask the accurate diffusion coefficient. Even in the case 
that the system has large local MFG, it will not be 
problematic if the root mean squared displacement by 
diffusion (vDt) is much larger than the periodicity of the 
local MFG, because the local MFG is expected to be av- 
eraged out during diffusion. On the contrary, in the case 
of systems with small diffusion coefficient, the local MFG 
can be remained and spin echo trains of magnetic reso- 
nance can be further attenuated by the local MFG as 
well as by the external MFG. 

So far, there are several kinds of possible approaches to 
cancel the effect of local MFG, that is to set up the exper- 
imental condition of G e 3> G; by production of large ex- 
ternal MFG by an anti-Hclmholtz superconducting mag- 
net or a fringe field of a superconducting magnet. By 
using anti-Hclmholtz superconducting magnet, one can 
obtain a MFG up to 200 T-m" 16 . On the other hands, 
by using superconducting fringe field, one can obtain up 
to 60 T-m _1 £. There are the other approaches to can- 
cel local MFG by using bipolar gradients£*£*2ii2i or rotary 
echosii. These methods are based on the idea of cance- 
lation of unwanted local MFG. On the contrary, we are 
interested in local MFG, which is thought to be char- 
acteristic to each systems particularly to solids, and it 
should be closely related to the molecular and electronic 
structure. 

There are a number of literatures available to calcu- 
late spin echos in restricted geometriesiiiiiSiiSiii. Most 
of these researches are devoted to describe pulse gradi- 
ent spin echo(PGSE) experiments. The spin echo ex- 
periment under a static field gradient is a kind of niche 
applications, and its theoretical treatment, namely, sol- 
vation of Bloch-Torrey equation^, is still lacking. Up 
to now, Axelrod and Sen have shown a solution of one 
dimensional Bloch-Torrey equation under the linear gra- 
dient condition by using an eigen mode technique^. In 
this Letter, we shall extend their treatment to spin echo 
experiments with two-dimensional lattice diffusion un- 
der arbitrary magnetic field gradient shape. Although 
the general theory of spectral density function for two- 
dimensional lattice diffusion is developed 19 , effects of 2D 
diffusion on nuclear spin echo experiments have not been 
investigated so far. 

The method developed here will be useful to inves- 
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tigate transport phenomena in low dimensional solids, 
such as movement of lattice defects (including soliton 
in disordered solids), dynamics of elementary excita- 
tion(including exciton and polaron), transport of charge 
carriers, and so on. Further, recent developments on 
magnetic resonance force microscopy (MRFM) show that 
one can perform in situ imaging with a scale of several 
tens of nanometers, where ultra high magnetic field gra- 
dients of over 10 3 ~ 10 4 T/m are generally used22i2i. 
Up to now, there are no publications available to obtain 
spin-echo measurements under a condition of such a huge 
MFG, which can be comparable to local magnetic field 
gradient with atomic scale. The theory presented in this 
Letter will be useful for analyses of such experiments as 
well. 



II. THEORETICAL 

The method presented here is based on the numeri- 
cal solution of the Bloch-Torrey equation by the Fourier- 
spectrum method^. Axelrod and Sen have proposed the 
methodology for calculating amplitudes of spin echos for 
a nucleus diffusing under a linearly inhomogeneous mag- 
netic field condition^. Although they have shown the 
solution with a linear magnetic field gradient by an eigen 
mode expansion method, the details of the numerical pro- 
cedure was not shown. In this article, we shall demon- 
strate a mathematical formulation of solvation of the 
Bloch-Torrey equation by the Fourier-spectrum method 
and explore spin echo amplitudes under various magnetic 
field gradients. 

So far, we have attempted to solve the differential equa- 
tion of diffusion equation in the real space by Crank- 
Nicolson method^ 3 -. However, we have suffered from the 
underflow of variables due to iterative multiplication in 
many times, of the time evolution operator for the spin 
echo experient. Furthermore, while the Crank-Nicolson 
method preserves the absolute stability against the size 
of each step in the real space, we observed the strong 
dependence of spin echo amplitudes on magnitude of in- 
terval of time step. 

Thus, these problems lead us to motivation for numer- 
ical simulation in the Fourier space concerning coordi- 
nates in real space. 

Below we shall describe the mathematical formulation 
for Bloch-Torrey equation with magnetic field gradient 
of an arbitrary wave form. 



A. Solution for ID Bloch-Torrey Equation with 
arbitrary magnetic field gradient 



The Bloch-Torrey equation is shown in the following: 

ijf(x)M, (1) 
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where the tilde denotes dimensionless parameter in or- 
der the problem to be treated as general. The effects 
of longitudinal relaxation (Ti) and intrinsic transverse 
relaxation(T2) in original Bloch-Torrey equation was ne- 
glected in the calculation. But the effect must of course 
be taken into account in practical applications. Dimen- 
sionless diffusion coefficient (Do), dimensionless gyromag- 
netic ratio(7), and one-dimensional coordinate(i) are de- 
fined as 
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(3) 
(4) 



f(x) is a functional form of ID MFG. 

Let us think the expansion of the Eq. (JJ) by a series of 
orthogonal functions, ^ k (x). 



dM 
~0T J 



DoT 



d 2 M 



ijf[f(x)M] 



(5) 



Here, we used linearity of orthogonal function transfor- 
mation. The second term of the right hand side of Eq.JSJ) 
is obtained by convolution: 



(6) 



Thus, one can write down the practical form of Eq.Q 
like the following: 



dM k (t) 
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-{nk) 2 D M k (t) - <7 ^ V k - p M p {t), (7) 

p=0 



Here. 



nl/2 

M fc (i) = F[M} = / M(x,t)V k (x)dx. (8) 

J -1/2 

Vk is a k-th coefficient for orthogonal function trans- 
form(for example, Fourier transform) of x: 



V k = 



1/2 



1/2 



f(x)V k (x)dx. 



(9) 



Further, the first term of the right hand side of Eq.Q) 
can be obtained by using the fact that second derivative 
in real space corresponds to the multiplication with k 2 . 



Here. 
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Let us think a basis set with taking account of the 
boundary condition of x: [-1/2:1/2]. Shown below is the 
eigen modes using sinusoidal(only cosine components) or- 
thogonal functions for the Laplacian, Ao = —V 2 , of the 
one dimensional diffusion equation, 



$ k (x) = C* fc cos(7rfc(i + -)), 



where 



C k 



1 if k = 
y/2 if k ^ 



(13) 



(14) 



and corresponding eigen values are = (irk) 2 . 

Since it is unnecessary to think the time ordering dur- 
ing the echo time, r, in transverse relaxation time mea- 
surements, one can take r as one step. 

The time evolution operator for the echo time, r, can 
be defined as 



U 4 



(W-ijB)r 



(15) 



Thus, the propagator for the two pulse Hahn echo^* can 
be express as 



u 2r = u-u+ = (u+yu^ 



(16) 



Similarly, that for the second echo of Carr-Purcell 
Meiboom-Gill(CPMG)SSi2& can be expressed as 



U 4t = U+U-U-U+, 
and that for n-th echo of CPMG, 



U 2 



TT-mod(n,2) j-i-floor(n,2) 
U 2t ' U 4t 



(17) 



(18) 



Therefore, The k-th Fourier component of the magneti- 
zation at a time, t = 2nr, can be calculated by 



M k (t) = U 2nT Mk(0). 



(19) 



Here, the initial magnetization, Mfc(O), can be obtained 
by using the initial condition of 4>{x): 



OO OO 

Y.MkiO^kii) = ^M fe (0)cos(7rfc(£+-)) 



fe=0 



k=0 

<t>(x). 



(20) 



The Fourier component, Mfc(0), of the initial magneti- 
zation, can be calculated by inverse Fourier transform 
of 4>(x). For example, in the case of homogeneous spin 
density, namely, in the case of the initial condition of 
4>(x) — 1 (constant), all the components are zero except 
k=0 since the inverse Fourier transform is a ^-function. 
Substituting the obtained -Mfe(O) into Eq.JEJJ, the k- 
th magnetization, Mk(t), at a time, t(= 2m"), can be 
computed. Finally, one can calculate the magnetization, 
M(t) by inverse Fourier transform of Mfc(i) and the sum- 
mation of each components which correspond to each 
components in real space; 



1/2 



M(t) = ^"MW)]. 



(21) 



£=-1/2 



All the figures in this article are plots of the relaxation 



exponent, log 



M(0) 
M(t=2nr) 



/n, as a function of Dq 



B. Two-dimensional Diffusion 

We extended the ID Bloch-Torrey equation to the case 
of anisotropic 2D lattice diffusion. For anisotropic 2D dif- 
fusion, the Bloch-Torrey equation is expressed as follows. 

dM ~ d 2 M - d 2 M 

-gf = D*o-Q^+Dy -^- i7 f(x,y)M. (22) 



The solution can be written as 

M(t) = C7 2 „ T M(0), 

where 



(23) 
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exp 



(-D V 2 -i-yf(x,y))dt) (24) 



Now we perform Fourier transform using the series of 
two dimensional orthogonal functions, ^ mn (x, y), 

V mn (x,y) = C mn ^ m {x)^ n {y) 

= C mn cos(Trm(x+-))cos(irn(y+-)) 

(25) 

C mn is a normalization constant and m, n are eigen 
mode values. 



1 (m = n = 0) 

V2 (m = or n = 0) (26) 

2 (m/0 and n ^ 0) 



On 



Let us consider the space Fourier transform of Eq. 

T[M t ] = D f[V 2 M}-i^[f(x,y)M}. (27) 

Again, T[M] stands for the 2D Fourier transform of 
magnetization M, 



M mn (t) = T[M] 



M{x,y,t)^ mn {x,y)dxdy (28) 



Here, M represents the magnetization in the Fourier 
space. We begin by the procedure for calculating 
F{V 2 M] of Eq.lgZJ. 



T[V 2 M] = T[M^\ + T[Myy] 



(29) 
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For the second-order derivative of the y component, 
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Similarly, for the i component, F[M XX \ = W X M, 
where, 



/-(tt-0) 2 

-(^-i) 2 



\ 





-(tt ■ m) 2 J 



(31) 



Thus, 



F[M ££ ] + T[Myy] = W X M + MWy 



(32) 



When m — n, both W x and W y are same diagonal 
matrices each other, so we can rewrite, 

T\M xi \ + T[Myy\ = W X M + MW X . (33) 
Now we introduce a Laplacian superoperator W based 



on commutator superoperator: 



,27 



WM = W X M + MW X 

= W X M x E + E x MW X , (34) 

where E represents the unity operator. This leads to the 
matrix representation, 



W = W X ®E + E®W X 



(35) 



where W x represents the transposed matrix of W x . ® 
denotes direct product of matrices. 

In the practical calculation, we assume the square lat- 
tice and adopt the identical number of eigen modes for x 
and y dimensions(m = n), and E = E. 



Now we go on to show the calculation of the second 
term, T\f(x, y)M], in Ea. (|27l) . The Fourier transforma- 
tion of this term can be rewritten using a convolution, 

T[f{x,y)M] = T[f(x,y)}*T[M] 

= J212 Vm -p< n -i M p><» ( 3? ) 

p q 

where V m . n represents (m, n)-th component of Fourier 
space of f(x, y). 

V m ,n= // f{x,y)^ mn (x,y)dxdy (38) 



Thus, the convolution superoperator, B, between 
J 7 [f(x,y)] and T[M mn ] is defined as 



F\J(x, y)M] = T[f{x, y)] * F[M\ -» BM. 



(39) 



Here, it is worth to note that the computational effort 
of the convolution in Ea. H37[l is proportional to 0(m 2 x 
(30) 77 , 2 ) ; so this calculation becomes very huge. To avoid this 
difficulty, one can employ the transform method utilizing 
fast-Fourier Transform (FFT) in oder to diminish the 
computational time 2 ". The computational effort of the 
transform method is proportional to O (m log mxn log n) . 
In the transform method, one calculate the elaborated 
formula: 



(40) 



jf[V 2 M] = WM 



(36) 



By the virtue of fast Fourier Transform(FFT) algorithm, 
one can reduce the number of multiplication to calculate 
convolution. Unfortunately, we could not find the su- 
peroperator for the convolution, then, we implemented 
the direct calculation of convolution by using a paral- 
lel programming technique of the Message Passing Inter- 
face(MPI) standard. 



III. RESULTS AND DISCUSSION 

A. Comparison with the result of the second-order 
cumulant expansion 

Figffl shows a comparison of the numerical solu- 
tion of Bloch-Torrey equation by the Fourier spectrum 
method, with the results of the second-order cumulant 
expansion(so-called gaussian phase approximation). The 
calculation was performed for the cases of a linear(/(x) = 
x) and parabolic(/(a;) = Ax 2 — 0.5) magnetic field gradi- 
ents in one-dimensional. According to Axelrod and Sen, 
the second cumulant expansion become a good approxi- 
mation for both the short time and motionaly averaging 
regimeaii. For the linear case, we obtained the exactly 
same behavior on the result of Fourier spectral method 
as the result by Axelrod and Senii. As for the localized 
regime, the relaxation exponent indicates strong depen- 
dence on the functional shape of a magnet field gradient. 
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B. Effect of Dimensionality of Diffusion 

A sinusoidal function may be a good candidate as a 
first step, for the expression of internal magnetic field 
gradient in a crystal. For example, the Peierls poten- 
tial is often approximated by the sinusoidal. Fig[3] shows 
the decay exponentials for the 8th CPMG echos un- 
der the condition of the sinusoidal gradients, shown in 
the figure, as a function of dimensionless diffusion co- 
efficient, Do (Dq = D x q for one-dimensional cases and 

D = (D x0 + DyQ ) x / 2 for two-dimensional cases). We 
took into account the lowest eigen modes up to 16th and 
used the 7 value of 10. For the so-called localization 
regime, the lower attenuation of decay exponential was 
observed for the 2D diffusion case than that for the ID 
case. It is also found that the dimensionality of diffusion 
is insensitive both at the short time and at the motional 
averaging regimes. At the motional averaging regimes, 
although a slight increase of echo attenuation occurs on 
2D cases, the decay exponent as a function of Dq shows 
the same slope irrespective of dimensionality of diffusion. 

C. Effect of Anisotropic Diffusion in a 2D lattice 

Now let us look at a situation that a system of inter- 
est shows the anisotropic diffusion from internal and/or 
external reasons. For example, putting the sample in the 
external electric field corresponds to the latter case. Fig|3| 
shows the effect of anisotropy of diffusion on CPMG de- 
cay exponent. One observe a double maximum profile 
of relaxation exponent for diffusion with an anisotropy 
of over 200. Two peaks of the relaxation exponent have 
some interesting features. 

First, while the main peak around the Dq value of unity 
does not nearly move its position with respect to changes 
in anisotropy, the second peak moves to positions with 
larger Dq values proportional to the anisotropy. Sec- 
ond, the ratio of the relaxation exponent at two max- 
ima is almost constant with respect to changes in the 
anisotropy of diffusion, but the ratio and the shape of 
peak are strongly dependent to functional form of the 
field gradient. Third, the differences in dimensionality of 



diffusion show much larger effect on relaxation exponent 
than functional form of MFG. By comparing Fig[3Ji) and 
b) with c) and d), one can realize that the relaxation 
exponents for two-dimensional MFG arc much larger 
than those for one-dimensional for the case of isotropic 
difrasion^o/A/O = 1). 

From the above results, it would be possible to de- 
termine the functional form of MFG by transverse relax- 
ation time measurements at localization regime when the 
anisotropy of diffusion is large. If the MFG is originated 
from internal one, for example, crystalline periodicity, 
Peierls potentialSiSLSi, Friedel oscillation^, and so on, 
transverse relaxation time measurements will become a 
novel tool for determination of the functional form of 
these internal properties. 



IV. CONCLUSION 

We have explored diffusion and internal local mag- 
netic field gradients in solids by nuclear magnetic res- 
onance(NMR) transverse relaxation time measurements. 
The Fourier-spectrum method was employed in order to 
solve the Bloch-Torrey equation with arbitrary magnetic 
field gradients in one- and two-dimensional lattice restric- 
tions. If there are no phase transitions or deformations 
of samples for a temperature region interested, variable 
temperature CPMG measurements with a fixed echo time 
have a possibility for determination of dimensionality of 
diffusion processes found in the samples. The method 
will be applicable to several research areas such as trans- 
port phenomena of charge carriers in ion conducting ma- 
terials, fluid dynamics in low dimensional porus materi- 
als, and dynamics of low dimensional elementary excita- 
tions in solids. 
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FIG. 1: Comparison of numerical solution of ID Bloch-Torrey 
equation by the Fourier spectrum(FS) method, with the re- 
sults of the second-order cumulant expansion(gaussain phase 
approximation;GPA). The calculation was performed for the 
cases of linear(/(a:) = x) and parabolic(/(:r) = 4a; 2 — 0.5) 
magnetic field grandients. The lowest eigen modes up to 64th 
were taken into account. 
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FIG. 2: Effect of dimensionality of diffusion on decay expo- 
nent of CPMG spin echo experiments. The eigen modes up 
to 16th is taken into account and the value of 7 is set at 10. 




FIG. 3: Effect of anisotropic diffusion on decay expo- 
nent of CPMG spin echo experiments. a) f(x,y) = 
0.5 cos(27ra) cos(27rj/), b) f(x,y) — 0.5 — cos(-7r:r) cos(-7rj/), c) 
f(x, y) = —0.5 cos(27ra;), and d) f(x, y) = x. 



